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Abstract 

General  delay  dynamical  systems  in  which  uncertainty  is  present  in  the  form  of  probability  measure 
dependent  dynamics  are  considered.  Several  motivating  examples  arising  in  biology  are  discussed.  A 
functional  analytic  framework  for  investigating  well-posedness  (existence,  uniqueness  and  continuous 
dependence  of  solutions),  inverse  problems,  sensitivity  analysis  and  approximations  of  the  measures  for 
computational  purposes  is  surveyed. 

Keywords:  Inverse  dynamics  problem,  Probability  distribution  function,  Sensitivity  analysis,  Biomedical 
systems 


1  INTRODUCTION 

The  purpose  of  this  presentation  is  to  survey  recent  as  well  as  forthcoming  results  in  our  research  efforts  on 
models  with  delays  and  hysteresis  where  probabilistic  uncertainty  is  present  in  a  significant  way.  While  we 
focus  our  motivation  here  on  examples  arising  in  biological  applications  (Banks  and  Bihari,  2001;  Banks  and 
Holte,  2003;  Banks  and  Potter,  2003;  Banks  and  Bortz,  2005a;  Banks  and  Bortz,  2005 6;  Banks  and  Davis,  to 
appear;  Banks  and  Pinter,  2005;  Banks  and  Allnutt,  to  appear;  Banks  and  Nguyen,  to  appear;  Banks  and 
Nguyen,  work  in  progress),  similar  systems  arise  in  other  applications  as  diverse  as  materials  (Banks  and 
Medhin,  submitted;  Banks  and  Webb,  1997a;  Banks  and  Webb,  1997&;  Banks  and  Pinter,  2004;  Banks  and 
Pinter,  to  appear;  Banks  and  Pinter,  submitted;  Banks  and  Pinter,  2005),  electromagnetics  (Banks  and 
Gibson,  2005;  Banks  and  Gibson,  to  appear),  physics,  communication  networks,  etc.  As  is  explained  here, 
there  are  a  wide  class  of  models  related  to  cellular  level  population  dynamics  that  lead  to  systems  of  the 
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x(t)  =  /  x(t  +  6)dP(d)  +  (1) 

J —r 

where  P  is  a  generally  unknown  probability  measure  that  must  be  estimated  from  aggregate  or  population 
level  (as  opposed  to  individual  level)  observations  or  data.  The  probability  measure  P  (which  we  shall 
also  refer  to  as  a  probability  distribution  when  no  confusion  results)  may  be  discrete,  absolutely  continuous 
(continuous)  or  a  combination  of  both.  In  addition  to  the  obvious  inverse  problems,  there  are  fundamental 
questions  related  to  modeling  of  uncertainty,  well-posedness,  sensitivity,  estimation  and  approximation.  The 
primary  goal  of  this  note  is  to  outline  a  theoretical  and  computational  framework  to  treat  these  problems. 

2  EXAMPLE  FROM  CELLULAR  PATHWAYS:  HIV  INFEC¬ 
TION 


Figure  1:  HIV  infection  pathway  in  acutely  infected  cells. 


Our  first  example  is  typical  of  delay  systems  that  arise  in  biochemical  pathways  and  cellular  level  kinetics 
of  drug  metabolism  as  well  as  other  synthesis  models.  In  (Banks  and  Holte,  2003;  Banks  and  Bortz,  2005a) 
the  authors  study  a  model  for  progression  of  Human  Immunodeficiency  Virus  (HIV)  at  the  cellular  level. 
The  model  involves  compartments  T,  A,  C,  and  V  for  in  vitro  blood  level  counts  in  mice  of  target  (CD4+) 
cells,  acutely  infected  cells,  chronically  infected  cells  and  active  viral  particles,  respectively.  Free  virus  V 
infects  target  cells  T,  transforming  them  into  acutely  infected  cells  A  which  at  some  time  later  become 
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chronically  infected  cells  C .  The  basic  pathway  for  infection  and  production  of  virus  for  acutely  infected 
cells  is  schematically  depicted  in  Figure  1.  For  models  in  which  the  individual  kinetics  for  loss  of  envelope  and 
capsid,  integration,  transcription,  and  assembly  are  not  detailed,  it  is  necessary  (see  (Banks  and  Holte,  2003)) 
to  include  a  delay  t\  from  the  time  of  infection  of  a  target  cell  T  until  it  first  produces  free  virus  V.  There 
is  also  some  delay  t2  before  an  acutely  infected  cell  A  becomes  a  chronically  infected  cell  C. 

Here  we  outline  a  brief  derivation  from  first  principles  (with  assumptions  based  on  the  biology)  that 
supports  a  mathematical  model  in  which  the  delays  are  treated  as  probabilistically  distributed  across  the 
population  of  cells  found  in  a  typical  in  vitro  culture. 

First  consider  the  delay  between  initial  acute  infection  and  the  cell  becoming  what  is  termed  a  chronically 
infected  cell  characterized  by  differences  in  cell  dynamics  (see  (Banks  and  Holte,  2003)).  It  is  biologically 
unrealistic  (and  unacceptable  in  the  modeling  to  biologists)  to  expect  an  entire  population  of  cells  to  simul¬ 
taneously  change  infection  characteristics  precisely  72  (t2  >  0)  hours  after  initial  viral  infection.  Therefore, 
one  might  suppose  that  the  delay  between  initial  acute  infection  and  chronic  infection  varies  across  the  cell 
population  (thus  mathematically  characterizing  the  intercellular  variability)  according  to  a  probability  distri¬ 
bution  P2  (which  is  not  assumed  to  necessarily  possess  a  density  p2  -  it  could  have  point  masses).  Denote  by 
C(t ;  t)  the  subpopulation  consisting  of  chronically  infected  cells  that  either  maintained  their  acute  infection 
characteristics  for  r  time  units  or  are  the  progeny  of  those  same  cells.  In  other  words,  for  some  r  >  0,  there 
exists  a  subpopulation  C(t ;  r)  of  the  chronically  infected  cells  which  either  spent  r  hours  as  acutely  infected 
cells  (before  converting  to  chronically  infected  cells)  or  are  descendants  of  cells  that  spent  exactly  r  hours 
as  acutely  infected  cells.  Thus,  the  rate  of  change  in  this  subpopulation  of  cells  is  governed  by 

C(t;  t)  =  (r„  -  5c  -  5X(t))C(t ;  r)  +  jA(t  -  r), 

where 

X(t)  =  A(t )  +  C{t)  +  T(t) 

is  the  total  number  of  CD 4+  cells  (infected  and  uninfected).  The  expected  value  of  the  population  of  chronic 
cells  is  given  by  integrating  with  respect  to  the  distribution  P2,  over  all  possible  delay  values,  obtaining 

pOO 

C(t)=£2[C(t,r)]  =  /  C(t-,T)dP2(T).  (2) 

Jo 

Here  the  parameters  rv,6c,  5  and  7  are  appropriate  rate  parameters  (for  details,  see  (Banks  and  Holte,  2003)). 
Therefore,  the  rate  of  change  in  the  total  population  of  chronic  cells  is  governed  by 

pOO 

C(t)  =  £2[C(t-,T)]  =  (rv-5c-5X(t))£2[C(f,T)]  +  j  A(t  -  t)cIP2(t) 

Jo 

C(0)  =  Co, 

where  Co  is  the  initial  condition  for  the  total  chronically  infected  cell  population. 

Next  consider  the  delay  between  viral  infection  and  viral  production  for  the  acutely  infected  cells  A{i). 
Again,  it  is  unreasonable  to  expect  the  entire  population  of  acutely  infected  cells  to  simultaneously  commence 
viral  production  t±  (t\  >  0)  hours  after  infection.  Suppose  that  the  delay  between  infection  and  production 
(for  acutely  infected  cells  A(t))  varies  across  the  population  with  probability  distribution  Pi  (again  we  do  not 
assume  absolute  continuity  of  the  associated  measure) .  We  also  partition  the  expected  total  viral  population 
V  into  those  virions  Va  produced  by  acutely  infected  cells  and  those  virions  Vq  produced  by  chronically 
infected  cells  so  that 

V  =  VA  +  Vc. 

Then  we  denote  by  Va  (i;  r)  the  subpopulation  of  virus  which  are  produced  by  an  acutely  infected  cell  r 
hours  after  being  infected.  Thus,  the  rate  of  change  in  this  subgroup  of  virions  is  governed  by 

VA{t\  t)  =  - cVA{t ;  r)  +  nAA{t  -  r)  -  niVA(t;  r)T(t). 

To  obtain  the  (expected)  number  of  virus  at  time  t  that  have  been  produced  by  acutely  infected  cells,  we 
must  integrate  with  respect  to  the  distribution  P\:  over  all  possible  delays 

pOO 

VA(t)  =  £i[VA{t’,T)]  =  /  VA(t-,T)dPi(r), 

Jo 
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which  yields  the  governing  equation  for  this  larger  subpopulation  of  virions 
VA{t)  =  £x[VA{t\T)\ 

poo 

=  -cVA(t)  +  nA  VA(t;  r)dPi(r)  -  nj[VA(t)T(t). 

Jo 

To  account  for  the  chronically  infected  cells  as  a  source  of  virions,  we  denote  by  Vq  the  subpopulation  of 
virions  produced  by  chronically  infected  cells.  Thus  the  equation  describing  the  rate  of  change  in  the  size  of 
this  subpopulation  is 

Vc(t)  =  -cVc(t)  +  ncC(t)  -  niVc{t)T(t ), 

where  the  expected  value  C  of  the  total  population  of  chronically  infected  cells  is  defined  in  equation  (2). 
Therefore,  the  governing  equations  for  the  total  population  of  virus  are  described  by 


V(t) 


no) 


Ex  [VA(t-T)  +  Vc{t)] 

- c{VA(t )  +  Vc{t))  -  ni(VA(t)  +  Vc(t))T(t)  +  ncC(t )  +  nA 

POO 

—cV(t)+nA  /  A(t  —  T)dPx(r)  +  ncC(t)  —  niV(t)T(t) 

J  o 

Vo, 


A(t  —  r)dPi(r) 


where  Vo  is  the  initial  condition  for  the  total  virions  population. 

Moreover,  we  assume  that  the  A  and  T  subclasses  have  no  subpopulation  structures,  and  are  therefore 
governed  by 

pOO 

A(t)  =  {rv- 8A- 8X{t))A{t) +nIV{t)T{t)  -  7/  A(t  -  r)dP2(T) 

Jo 

^(0)  =  A) 

T(t)  =  (ru-6u-SX(t)-niV(t))T{t)  +  S 
T(  0)  =  T0, 


with  initial  conditions  Aq  and  To-  Note  that  in  equation  (3),  the  rate  term  with  the  delay  (representing  the 
delayed  conversion  of  A  to  C)  is  simply  the  negative  of  the  corresponding  delay  rate  term  in  equation  (3). 

Finally,  we  make  the  change  of  variables  Pj(£)  =  Pj(— so  that  the  distributions  are  now  defined  on 
(—oo,0)  instead  of  (0,oo)  (we  do  this  to  be  consistent  with  the  standard  notation  in  the  FDE  literature), 
and  obtain  the  system 


m 

A(t) 

C{t) 

m 


—cV(t)  +  nA  f  A(t  +  d)dPx(O)  +  ncC(t)  -  mV(t)T(t) 

J  —OO 


(rv  -  SA  -  6X(t))A(t)  +  niV(t)T(t)  —  7  f  A(t  +  0)dP2(6l) 


(rv  -SA-  6X(t))C(t)  +  7  /  A{t  +  0)dP2{0) 

J  —  OO 

(r„  -  8a  -  6X(t)  -  niV(t))T(t )  +  S, 


(3) 


which  is  a  vector  system  of  the  form  (1).  Special  cases  of  such  systems  include  those  in  which  the  probability 
measures  are  defined  on  some  finite  interval  Q  =  [— r,  0]  of  possible  delay  values  6  as  in  (1). 

This  model  was  successfully  used  (see  (Banks  and  Holte,  2003))  to  describe  in  vitro  mice  data  from  Dr. 
Michael  Emerman’s  lab  at  Fred  Hutchinson  Cancer  Research  Center.  Using  inverse  problem  methodology 
and  statistical  analysis  it  was  shown  that  improvement  of  fit  to  data  by  inclusion  of  the  delays  Ti  or  P\ 
is  statistically  significant  while  inclusion  of  delays  r2  is  less  important.  Indeed,  it  was  found  that  the 
experimental  data  could  not  be  properly  fit  with  an  ODE  version  (i.e. ,  with  the  delays  omitted)  of  the 
model. 
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3  EXAMPLE  FROM  A  VACCINE  PRODUCTION  MODEL 


A  second  class  of  models  (Banks  and  Allnutt,  to  appear)  that  illustrate  the  type  of  problems  focused  on  here 
involve  the  use  of  shrimp  grown  in  production  “raceways”  (essentially  large  growth  chambers  where  envi¬ 
ronmental  factors  such  as  temperature,  oxygen,  nutrient  levels,  etc.,  can  be  carefully  controlled)  artificially 
infected  to  efficiently  produce  large  quantities  of  an  associated  vaccine.  Scientifically,  this  entails  recruiting 
the  biochemical  machinery  in  an  existing  biomass  for  the  production  of  a  vaccine  or  antibody  by  infection 
using  a  virus  carrying  a  passenger  gene  for  the  desired  antibody  response. 

While  the  model  of  (Banks  and  Allnutt,  to  appear)  is  specific  to  virus  growth  and  vaccine  production  in 
shrimp,  the  implications  for  other  crustaceans  are  obvious.  And  of  course  the  shrimp  models  we  investigate 
can  serve  as  a  foundation  for  understanding  viral  progression  in  other  species  important  to  marine  agriculture. 
The  mathematical  goal  is  to  model  a  system  wherein  one  uses  shrimp  as  a  scaffold  organism  to  produce 
biological  countermeasures.  In  such  a  system  one  might  first  stock  shrimp  postlarvae  and  allow  them  to 
grow  normally  in  the  controlled  environment.  Then  one  infects  them  with  a  recombinant  viral  vector  (e.g., 
recombinant  Taura  Syndrome  virus  or  rTSV  in  the  example  developed  in  (Banks  and  Allnutt,  to  appear)) 
expressing  a  foreign  antigen,  resulting  in  vaccine  production  in  live  infected  shrimp. 

To  mathematically  demonstrate  the  feasibility  of  this  approach  one  considers  a  hybrid  model  of  the 
shrimp  biomass/countermeasure  production  system  which  has  two  components:  biomass  production,  and 
production  of  countermeasure  (antibody/vaccine).  The  output  of  the  biomass  production  model  is  input  to 
the  vaccine  production  model.  For  initial  investigations  the  amount  of  vaccine  produced  is  assumed  equal 
to  the  total  infected  biomass.  Thus,  the  vaccine  production  model  will  essentially  follow  the  course  of  the 
viral  dynamics  in  the  shrimp. 

The  effort  requires  modeling  the  dynamics  of  shrimp  at  the  population  level.  In  such  models  ignoring 
structure  in  constructing  mathematical  models  for  the  dynamics  of  shrimp  is  unrealistic,  since  shrimp  have 
size  dependent  characteristics  as  well  as  responses  to  external  environment.  An  appropriate  beginning 
model  is  based  on  the  classical  McKendrick-von-Foerster/Sinko-Streifer  size-structured  population  equations 
(Kot,  2001;  Metz  and  Diekmann,  1986)  with  mass  as  the  structure  variable,  i.e.,  one  equates  the  size  variable 
with  the  mass  in  the  model. 

While  there  appears  to  be  a  dearth  of  literature  on  modeling  epidemics  in  shrimp  populations,  in  (Lotz  and 
Breland,  2003)  the  authors  develop  a  non-structured  five  compartment  epidemic  model  of  TSV  that  includes  a 
Reed-Frost  transmission  process  in  closed  populations  of  shrimp  ( Litopenaeus  vannamei).  However,  structure 
can  play  an  important  role  in  the  study  of  viral  epidemiology  in  shrimp.  Moreover,  experimental  results 
(Hasson  and  White,  1999)  suggest  that  the  mortality  rate  in  acutely  infected  shrimp  depends  on  the  length  of 
time  that  the  shrimp  remain  acute.  Also,  individuals  in  the  latent  phase  have  varying  residency  times  before 
they  progress  into  the  acute  phase.  To  incorporate  all  of  these  features,  the  authors  of  (Banks  and  Allnutt,  to 
appear)  attempted  to  model  the  progression  of  TSV  in  shrimp  in  a  system  of  delay  PDEs.  However,  it  is 
difficult  to  correctly  account  for  the  different  residency  periods  of  individual  shrimp  in  this  fashion  as  the 
size  of  the  shrimp  is  a  function  of  time.  Instead  of  tracing  back  in  time  to  incorporate  delays,  a  different 
approach  involves  recording  the  variable  residency  times  in  the  different  stages  by  introducing  a  new  variable 
which  one  calls  the  class  age  of  an  individual.  The  class  age  of  an  individual  in  a  given  stage  represents  the 
length  of  time  that  the  individual  spends  in  that  stage  and  serves  as  a  surrogate  for  time  delays.  A  similar 
approach  has  been  used  to  investigate  a  linear  cell  population  model.  In  such  models,  cells  are  assumed 
capable  of  simultaneous  proliferation  and  maturation  where  in  the  proliferating  phase,  cells  are  committed 
to  undergo  cell  division  some  time  units  after  entering  this  phase  (Adimy  and  Pujo-Menjouet,  2003). 

In  the  vaccine  production  stage  the  shrimp  are  infected  by  distributing  chopped  dead  shrimp  infected 
with  a  recombinant  virus  evenly  throughout  the  raceway.  This  transfected  biomass  is  sufficiently  large  so 
that  most  of  the  shrimp  can  be  infected  in  a  short  period  of  time,  such  as  one  day.  There  are  other  modes  of 
transmission  of  virus  in  shrimp,  such  as  cohabitation  with  infected  shrimp  that  may  be  shedding  the  virus  into 
the  surrounding  medium  (waterborne  infection).  However,  compared  to  the  probability  of  shrimp  becoming 
infected  via  ingestion,  these  modes  of  transmission  can  be  assumed  (reasonably  for  a  first  investigation)  to 
be  negligible.  Hence  one  might  only  assume  infection  via  ingestion  of  dead  transfected  biomass.  It  might 
be  further  assumed  that  all  the  shrimp  have  an  equal  chance  of  becoming  infected  by  eating  the  infected 
biomass.  A  reasonable  time  interval  for  infection  is  7  to  10  days.  From  (Lotz  and  Breland,  2003)  and  (Dhar 
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and  Walker,  2004)  one  knows  that  during  this  time  interval  almost  no  shrimp  progress  into  the  chronic  state. 
Therefore  it  is  reasonable  to  consider  the  following  three  compartment  states:  susceptible  (S),  latently 
infected  (L)  and  acutely  infected  (A)  in  a  model.  In  this  model,  it  is  assumed  that  shrimp  will  become 
instantly  infected  (i.e.,  progress  into  latent  state)  as  soon  as  they  ingest  some  of  the  infected  biomass.  As  we 
have  noted  earlier,  however,  experimental  observations  suggest  that  there  exists  a  temporal  delay  between 
the  initial  latent  infection  and  initial  acute  infection  (Hasson  and  White,  1999).  Moreover,  it  is  biologically 
unrealistic  to  expect  all  members  of  the  shrimp  population  to  progress  into  the  acute  phase  at  a  fixed  number 
of  days  after  initial  latent  infection.  In  addition  the  shrimp  in  the  acute  phase  have  varying  mortality  rates 
because  of  the  different  times  that  they  progress  into  the  acute  phase  and  also  due  to  the  differences  in 
genetic  make-up  of  the  host.  As  we  have  already  noted,  it  is  difficult  to  account  for  the  class  age  history 
(i.e.,  the  length  of  time  that  shrimp  spend  in  a  state)  of  shrimp  in  a  particular  (latent  or  acute)  state  using 
a  system  of  delay  PDE’s  with  only  size  as  the  structure  variable.  This  is  because  it  is  not  obvious  how  to 
correctly  represent  the  integral  involving  the  delay.  As  an  alternative,  in  order  to  model  variable  residency 
times  one  may  keep  track  of  the  class  age  and  the  size  of  shrimp  by  incorporating  both  size  structure  and 
class  age  structure  into  the  latent  and  acute  states. 

Based  on  experimental  findings,  it  is  reasonable  to  assume  that  there  is  a  positive  probability  that  shrimp 
can  stay  in  each  the  latent  and  acute  state  for  more  than  7  to  10  days.  Thus  one  can  assume  that  the  class 
age  interval  for  both  states  is  the  same  as  the  time  interval  Tv  that  we  consider  in  our  model.  Note  that  all 
shrimp  from  the  biomass  production  raceway  are  healthy;  there  are  no  latently  infected  or  acutely  infected 
shrimp  in  the  raceway  at  time  t  =  0.  We  also  know  that  shrimp  in  the  acute  state  stop  growing,  which 
means  that  the  growth  rate  in  this  state  is  g  =  0. 

Based  on  the  above  discussions,  the  vaccine  production  model  developed  in  (Banks  and  Allnutt,  to 
appear)  is  given  by 

dtS(x,t)  +  dx(gs  (x)S(x,t))  +  ms  (x)S(x,t)  =  —A  S(x,t), 

dtL(x,  t,  9)  +  dx(gL(x)L(x ,  t,  9))  +  dgL(x,t ,  9)  +  mL(x)L(x1t ,  9)  =  — pL{9)L(x,t1 9 ), 

dtA(x,  t,  9)  +  dgA(x,  t,  9)  +  mA(9)A( x,  t,  9)  =  0, 

L(x,  t,  0)  =  XS(x,  t),  ^ 

A(x,t,  0)  =  [  L(x,t,9)dPL(9), 

Jo 

S{0,t)=0,  L(0,t,  9)  =  0,  A(0,t,9)=0, 

S(x,  0)  =  S°(x),  L{x ,  0, 9)  =  0,  A(x,  0, 9)  =  0, 

where  (x,  t,  9)  £  [xmjn,  .Tmax]  x  [0,  Ty\  x  [0,  Ty\.  In  the  above  S(x,  t)  denotes  the  density  of  individuals  having 

d 

mass  x  at  time  t  and  <9t  =  — .  The  function  L(x,t,9)  denotes  the  density  of  individuals  having  mass  x 

at  time  t  that  have  spent  9  days  in  the  latent  state  ,  whereas  the  function  A(x,  t,  9)  denotes  the  density  of 
individuals  having  mass  x  at  time  t  that  have  spent  9  days  in  the  acute  state.  The  quantity  gs( x)  denotes 
the  growth  rate  of  individuals  in  the  susceptible  state  ,  gL(x )  denotes  the  growth  rate  of  individuals  in  the 
latent  state.  The  function  ms (x)  denotes  the  mortality  rate  of  individuals  in  the  susceptible  state,  and  the 
function  mL(x)  denotes  the  mortality  rate  of  individuals  in  the  latent  state,  and  mA(6)  denotes  the  mortality 
rate  of  the  shrimp  that  spend  9  days  in  the  acute  state.  The  latent  to  acute  probability  density  rate  function 
Pl(9)  =  T'l(0)  defined  for  9  £  [0,7V]  denotes  the  rate  at  which  the  shrimp  in  the  latent  state  that  have 
spent  9  days  in  the  latent  state  become  acutely  infected,  while  the  quantity  A  denotes  the  infection  rate  due 
to  ingestion  of  chopped  infected  shrimp.  Finally  S°(x)  denotes  the  initial  population  density  of  susceptible 
shrimp  produced  from  the  biomass  production  model. 

We  note  that  this  is  again  a  probability  distribution  (Pl)  dependent  dynamical  system  (in  this  case  a 
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complicated  system  of  partial  differential  equations)  for  which  the  distribution  Pl  must  be  estimated  in  some 
type  of  inverse  problem. 


4  PROBLEM  FORMULATION  FOR  DISTRIBUTION  DEPEN¬ 
DENT  DYNAMICS 


In  both  the  examples  cited  above  as  well  as  in  many  others,  a  major  effort  involves  estimation  of  the 
probability  distributions  P  from  data.  A  typical  inverse  problem  consists  of  minimizing  the  output  least 
squares  criterion 


J(P)  =  £  | x(ti\  P )  -  di 

2  =  1 


(5) 


where  {di}  is  the  data  and  x(t;  P)  represents  the  solution  to  the  distribution  dependent  dynamics  such  as 
(1),  (3),  or  (4).  The  minimization  is  to  be  carried  out  over  the  space  V{Q)  of  probability  measures  defined 
on  a  set  Q  of  possible  parameter  values.  For  example,  in  the  systems  presented  above,  the  delay  times  or 
residency  times  are  restricted  to  some  finite  interval  Q  =  [ — r-,  0]  or  Q  =  [0,Ty],  respectively.  In  addition 
to  such  inverse  or  estimation  problems,  we  are  also  concerned  with  other  questions  for  problems  that  have 
distribution  dependent  dynamics  including  the  existence  and  uniqueness  of  solutions  to  the  dynamical  system, 
continuous  dependence  of  solutions,  sensitivity  of  solutions  with  respect  to  the  probability  distributions,  and 
numerical  approximations.  In  order  to  deal  with  such  questions  (either  theoretical  or  computational)  one 
needs  a  topology  on  the  measure  space  V(Q).  Indeed,  one  needs  a  number  of  items  to  develop  theoretical 
and  computational  foundations  including 


(i)  A  topology  on  V(Q), 

(ii)  Continuity  of  P  — >  J(P)  in  this  topology, 

(iii)  Compatible  compactness  results  (for  well-posedness), 

(iv)  Approximations  in  this  topology  (for  computations). 


It  is  fortunate  that  probability  theory  offers  significant  conceptual  help  toward  a  possible  complete,  tractable 
computational  methodology  (Billingsley,  1968).  The  primary  tool  is  the  Prohorov  metric,  which  can  be 
formally  defined  as  follows.  Suppose  (Q,d)  is  a  complete  metric  space.  For  any  closed  subset  F  C  Q  and 
e  >  0,  define 

Fe  =  {q  G  Q  :  d{q,  q)  <e,q€  F}. 

Define  the  Prohorov  metric  p  :  V(Q)  x  V(Q)  — >  R+  by 


p(Pi,P2)  =  inf {e  >0:P1[F}<  P2[Pe }  +  e,  F  closed,  F  C  Q}. 


This  can  be  shown  to  be  a  metric  on  V(Q)  that  satisfies 

(a)  (V(Q),p)  is  a  complete  metric  space, 

(b)  If  Q  is  compact,  then  (P(Q),p)  is  a  compact  metric  space. 

We  note  that  the  definition  of  p  is  not  intuitive.  It  is  not  clear,  for  example,  what  “P*,  — >  P”  in  p  metric 
means.  One  finds  the  following  important  characterization  (Billingsley,  1968). 

Theorem  1  Given  Pk,  P  £  V(Q),  the  following  convergence  statements  are  equivalent: 

(i)  p(Pk,P)  —>  0, 

(ii)  /  fdPk{q)  — >  /  fdP(q)  for  all  bounded,  uniformly  continuous  f  :  Q  — >  R1, 

JQ  JQ 
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(iii)  Pk[A]  — >  P[A\  for  all  Borel  sets  A  c  Q  with  P[dA]  =  0. 

Additional  useful  results  include: 

•  Let  C*b(Q)  denote  the  topological  dual  of  CB(Q ),  where  CB(Q )  is  the  usual  space  of  bounded  continuous 
functions  on  Q  with  the  supremum  norm.  If  we  view  V{Q)  C  CB(Q),  convergence  in  the  p  topology  is 
equivalent  to  weak*  convergence  in  V(Q)7 

•  Convergence  in  the  p  metric  is  equivalent  to  convergence  in  distribution. 

Considering  (ii)  of  Theorem  1,  it  is  readily  argued  that  the  dynamics  of  systems  such  as  (1)  (and  (3), 
(4))  are  p-continuous  in  P  on  V{Q).  Standard  arguments  from  the  theory  of  differential  equations  can  then 
be  used  to  argue  that  the  mapping  P  — >  x(t;  P)  is  also  continuous  on  {V{Q),p).  This  yields  immediately 
that 

n  2 

p^j(p)  =  '52\z(ti-,p)-di 

i—1 

is  continuous  in  the  p  topology.  Continuity  of  P  — >  J(P)  and  compactness  of  V(Q)  (each  with  respect  to 
the  p  metric)  allows  one  to  assert  the  existence  of  a  solution  to  min  J(P)  over  P  £  V(Q). 

As  we  shall  see  below,  the  Prohorov  metric  is  also  fundamental  to  development  of  a  sensitivity  theory  as 
well  as  “finite  element”  type  approximation  schemes  for  the  systems  of  interest  here. 


4.1  Abstract  Formulation 

Before  we  present  some  theoretical  results  on  a  class  of  distribution  dependent  problems,  we  illustrate  how 
to  derive  an  abstract  Cauchy  formulation  in  a  complex  Banach  space  for  a  typical  delay  system  example  such 
the  HIV  model  above.  First,  we  assume  the  HIV  system  derived  in  Section  2  only  depends  on  absolutely 
continuous  (continuous)  distributions;  that  is,  dPi{r )  =  pj(r)dT  for  i  =  1,2  where  Pi(r)  are  probability 
densities.  A  more  general  framework  for  discrete  distributions  and  mixed  (a  combination  of  continuous  and 
discrete)  distributions  is  discussed  in  the  next  section  where  one  assumes  the  form 

k 

-P(r)  =  (6) 

i=  1 


and 


k 

p(T )  =  Ep^t)  + 

i=  1 


(7) 


Here  Ar  is  the  Dirac  measure  with  atom  (mass)  at  r  for  the  discrete  and  mixed  measures,  respectively. 

Returning  to  the  HIV  model  we  let  v  =  (V,  A,  C,  T)t  and  x(t)  =  (v(t),vt)  £  X  =  R4  x  L2(— r,  0;K4). 
For  0  <  r  <  oo,  we  denote  the  parameter  space  M  =  L2(— r,  0)  x  L2(— r,  0)  and  Mc  =  {{pi,P2)  €  M\p\, 

P2  >  0  and  /  pi(r)dT  =  1}.  Then  the  HIV  system  derived  in  equation  (3)  can  be  rewritten  as  an  abstract 

J —r 

Cauchy  problem 


x(t)  =  Ax(t)  +  /2(t),  t>  0 

z(0)  =  x0,  (8) 

where  fi (t)  =  ((0, 0, 0,  S)T ,  0)  £  X ,  and  xo  =  (p,  (j>)  £  X.  Here  A  is  a  nonlinear  operator  such  that  A  : 
V{A)  Cl-tl  and  A(r),  <j>)  =  (^L(rj,  (f)  +  /i(p),  where  V(A)  =  {(rj  ,</>)£  X  \  <f>  £  H1(—r,  0;  R4)  and  p  = 
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0(0)}.  Furthermore,  for  (r;,0)  el4x  L2(— r,  0;  R4), 

L(v,  0)  = 


—c  0  nc  0 

0  rv  —  5  a  0  0 

0  0  rv  —  Sc  0 

0  0  0  ru  —  Su 


+nA [<5(i,2)]  (4,4)  /  0(T)pi(r)dr 
J  —r 


+7(  [<5(3,2)]  (4,4)  -  [S(2,2)](4,4))  /  0(r)p2(T)dr, 


-0771^4 

4 

V2  +  amm 

i= 2 

/l(’,)  =  -set*)* 

i= 2 
4 

-<5(5^  Vi)7)*  -  Otr)\T]i 
i= 2 

where  [<5(*,j)](4,4)  is  a  4  x  4  matrix  with  a  one  in  the  (i,j)th  component  and  zeros  everywhere  else.  In  (Banks 
and  Bortz,  2005a),  (Banks  and  Holte,  2003)  the  mass  action  product  nonlinearities  in  fi  are  replaced  by 
saturating  nonlinear  functions  -  see  the  definition  of  f1  in  (Banks  and  Bortz,  2005a),  (Banks  and  Holte,  2003). 

Once  an  abstract  Cauchy  formulation  is  constructed,  existence  and  uniqueness  of  solutions  for  equation 
(8)  follow  from  the  results  in  (Banks  and  Bortz,  2005 &;  Banks  and  Nguyen,  to  appear)  which  are  summarized 
in  Section  4.2  below.  Moreover,  theoretical  results  in  (Banks  and  Nguyen,  to  appear)  also  provide  continuous 
dependence  of  solutions  along  with  the  derivation  of  the  sensitivity  function  for  general  nonlinear  ordinary 
differential  equations  (ODEs)  in  a  Banach  space.  Here  we  only  show  the  construction  of  the  abstract  Cauchy 
problem  for  delay  systems  with  continuous  probability  measures.  However,  an  abstract  Cauchy  formulation 
for  delay  systems  with  discrete  and  mixed  distributions  of  the  form  (6)  and  (7),  respectively,  can  also 
be  constructed  using  similar  concepts  but  with  different  parameter  spaces.  These  theoretical  results  and 
associated  parameter  spaces  of  probability  measures  are  the  focus  of  the  next  section. 

4.2  Theoretical  Results 

In  this  section  we  recall  theoretical  results  that  treat  delay  systems  with  absolutely  continuous  (continuous) 
distributions.  Interested  readers  can  find  more  details  on  the  theories  and  the  proofs  in  (Banks  and  Nguyen, 
to  appear).  As  one  can  see  from  the  previous  section,  time  delay  problems  with  absolutely  continuous 
(continuous)  distributions  are  a  special  case  of  an  abstract  nonlinear  ODE  where  the  state  space  is  a  general 
Banach  space  X  and  the  parameter  space  A4C  is  a  convex  subset  of  a  Banach  space  A4  such  as  A4  =  L2  x  L2. 
Therefore,  consider  a  general  nonlinear  ODE  of  the  form 

x{t)  =  (9) 

where  /  :  K+  x  X  x  A4  — >  X  and  X  and  A4  are  complex  Banach  spaces.  We  define  the  successive 
approximations  for  system  (9)  to  be  the  functions,  x  ,  a;1, ...,  given  recursively  by 

x°(t,tQ,x  o,/r)  =  xQ, 

xk+1(t,t0,x0 ,/z)  =  x0+  f(s,xk(s,t0,x0,p),iJ,)ds, 

Jt0 

for  k  =  0, 1,2, ....  Then  one  can  establish  the  following  theoretical  results. 
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Lemma  1  (Existence  and  Uniqueness  of  Solutions)  Let  f  :  R+  x  X  x  M  — >  X  be  continuous  and 

\f{t,xi,n)  -  f(t,x  2,n)\  <C\X!-  x2\ 

for  some  constant  C  >  0.  Then  the  successive  approximations  xk  converge  uniformly  for  t  G  [to,T]  to  a 
unique  solution  x  of  (9)  such  that  x(to,  to,  xo,  p)  =  xq- 

Lemma  2  ( Continuous  Dependence  of  Solutions  on  Parameters)  Let  f  G  C[M+  x  X  x  Ad,  X ]  and  for  p  =  po, 
let  x(t,  to,  xq,  po)  be  the  solution  of 


x  =  f(t,x,p0),  x{t0)=x  o, 

existing  on  [to,T].  Assume  further  that 


lim  f(t,x,p )  =  f{t,x,p0), 

*A»o 

uniformly  in  (t,  x)  and  for  (t,  x\ ,p),  (t,  x2l  p)  €  R+  x  A'  xM, 

\f(t,X!,p)  -  f(t,x2,p)\  <  C\xx  -  x2\ 


for  some  constant  C  >  0.  Then  the  differential  system 

x  =  f(t,x,p),  x(t0)  =  x  o, 

has  a  unique  solution  x(t,  to,  Xo,  p)  satisfying 

lim  x(t,t0,x0,p)  =  x(t,tQ,x0,p0),  te[t0,T}. 

Even  though  the  results  given  here  are  under  the  assumed  global  Lipschitz  condition,  similar  results  can  also 
be  established  under  the  weaker  assumptions  of  a  local  Lipschitz  condition  and  /  being  dominated  by  an 
affine  function.  We  let  B( X,  Y)  denote  the  space  of  bounded  linear  operators  from  X  onto  Y  and  summarize 
a  sensitivity  theory  for  delay  systems  with  absolutely  continuous  (continuous)  measures  in  Theorems  2  and 

3. 


Theorem  2  Suppose  the  function  f(t,x,p)  of  (9)  has  continuous  Frechet  derivatives  fx(t,x,p)  with  respect 
to  x  and  ffft,x,p)  with  respect  to  p  with 


I  fx{t,x,p)\  <  M0  and  {f^t,  x,  p)\  <  Mx 


for  some  constants  Mq  >  0  and  M\  >  0.  Then  the  Frechet  derivative  y(t) 
in  B{M,X)  satisfying  the  equation 


d_ 

dp 


x(t,to,xo,  p)  exists  with  y(t) 


y(t )  =  fx(t,x(t,t0,x0,p),p)y(t)  +  fu.(t,x(t,t0,x0,p),p), 

y(to)  =  o, 


for  t  >  t0. 

With  the  parameter  space  of  probability  density  functions  Adc,  which  is  a  convex  subset  of  a  Banach  space 
A4 ,  the  sensitivity  theory  (Theorem  2)  above  can  also  be  applied  using  directional  derivatives  instead  of  the 
Frechet  derivative.  However,  it  is  shown  in  (Banks  and  Nguyen,  to  appear)  that  the  directional  derivative 
of  a  continuous  function  g  is  the  Frechet  derivative  on  A4  restricted  to  q  —  p  where  p,q  G  M. c. 

In  order  to  accommodate  delay  problems  with  Dirac  (discrete)  measures  or  measures  with  a  continuous 
component  and  a  saltus  component,  the  theoretical  results  above  are  extended  to  a  general  convex  metric 
space.  This  is  necessary  because  the  parameter  space  associated  with  the  Prohorov  metric  is  no  longer 
a  Banach  space  but  only  special  case  of  a  convex  metric  space  (Afi,  dxj-  For  discrete  measure  delay 
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systems  where  the  measures  are  defined  in  (6),  the  parameter  space  {Mi,  dj^i)  is  normally  chosen  to  be 
a  topology  with  the  Prohorov  metric  {V{Q),p).  Although  the  Prohorov  metric  is  not  conceptually  easy  to 
use,  it  generates  a  similar  topology  to  the  weak  L2  topology  (e.g.,  see  (Banks  and  Pinter,  2005))  which  is 
of  course  the  same  as  the  weak*  topology  in  this  case.  Therefore,  the  Prohorov  metric  may  be  applied  in 
numerical  approximation  in  distribution  dependent  problems  taking  advantage  of  its  relation  in  convergence 
to  the  weak  L2  convergence.  For  delay  systems  with  mixed  measures  as  defined  in  (7),  the  parameter  space 
can  be  based  on  a  combination  of  the  Prohorov  metric  topology  and  the  weak  L 2  topology  for  compatibility. 
Therefore,  Banks,  Dediu  and  Nguyen  have  extended  the  theoretical  results  mentioned  above  to  the  case 
where  the  parameter  space  is  a  convex  metric  space.  Let  {Mi,  dMi)  denote  a  general  convex  metric  space 
with  distance  dj^h  and  X  denote  a  general  complex  Banach  space.  Consider  a  general  nonlinear  abstract 
ODE 


x(t)  =  f{t,x(t),n  i),  x(0)  =  0, 


(10) 


where  /  :  R+  x  X  x  Mi  —>  X  is  continuous  in  all  three  variables  and  Frechet  differentiable  in  x.  Here 
the  solution  x  £  X  and  the  parameter  pi  £  Mi-  The  conditions  for  and  statement  of  existence  and 
uniqueness  of  solutions  of  equation  (10)  along  with  continuous  dependence  of  solutions  for  the  general  convex 
metric  parameter  space  are  similar  to  those  for  the  situation  detailed  above  where  the  parameter  space  is 
a  general  complex  Banach  space;  therefore,  those  theoretical  results  are  not  repeated  here.  When  deriving 
the  sensitivity  theory  for  the  convex  metric  parameter  space  case,  the  directional  derivative  is  used  instead 
of  the  Frechet  derivative  with  respect  to  the  measures. 

Given  any  two  arbitrary  points  p\,  v  £  {Mi,  dj^ q),  we  define  the  directional  derivative  Sf{t,  x,  p£,  v  —  p i) 
of  /  at  pi  in  the  direction  v  —  pi  to  be  the  value  of  the  limit 


lim 

e — >0 
e>0 


f{t,  x,  m  +  e{is  -  pi))  -  f(t,  x,  m) 
e 


8f(t,x,ni\v-  m), 


provided  this  limit  exists  in  X .  A  sensitivity  theory  for  a  convex  metric  parameter  space  is  stated  next;  more 
details  and  proofs  can  be  found  in  (Banks  and  Nguyen,  work  in  progress). 

Theorem  3  Suppose  the  function  f{t,x,p i)  of  (10)  has  a  Frechet  derivative  fx{t,x,p i)  with  respect  to  x 
such  that  fx  £  C[R+  x  X  x  Mi,B{X,X)\  and  \fx(t,x,  pi)\  <  M0  for  some  constant  M0  >  0.  More¬ 
over,  assume  f  also  has  a  continuous  directional  derivative  8f(t,x,pi;v  —  fii)  with  respect  to  pi  in  the 
direction  of  {v  —  pi)  such  that  \8f(t,x,pi\v  —  /ii)|  <  Mi  where  Mi  >  0.  Then  the  directional  derivative 
y(t)  =  Sx{t,  pi]  v  —  pi)  exists,  with  y  :  R+  x  X  x  M i  — >  X ,  and  y  satisfies  the  equation 


y{t )  =  fx{t,x{t,t0,x0,pi),pi)y{t)  +  5f{t,x{t,t0,x0,pi),pi\ v  -  pi), 

y{to )  =  o.  (li) 


Having  presented  theoretical  results  to  deal  with  delay  differential  equations  where  the  time  delay  is 
distributed  with  different  types  of  probability  measures  (i.e,  absolutely  continuous,  continuous,  discrete  and 
mixed  measures),  we  next  discuss  some  numerical  approximation  issues  for  this  class  of  problems. 


4.3  Approximation  Issues 

We  first  note  that  even  when  the  parameter  set  Q  is  finite  dimensional,  the  metric  space  ( V{Q ),  p)  is  infinite¬ 
dimensional  and  hence  one  must  use  finite-dimensional  approximations  to  obtain  tractable  computational 
algorithms.  To  this  end,  one  may  prove  (see  (Banks  and  Bihari,  2001)) 

Theorem  4  Let  Q  be  a  complete,  separable  metric  space  with  metric  d,  B  the  class  of  all  Borel  subsets  of 
Q  and  V{Q)  the  space  of  probability  measures  on  {Q,B).  Let  Q o  =  be  a  countable,  dense  subset  of 

Q.  Then  the  set  of  P  £  V{Q)  such  that  P  has  finite  support  in  Qo  and  rational  masses  is  dense  in  P{Q)  in 
the  p  metric.  That  is, 

k  k 

Vo{Q)  =  {P  £  V{Q)  :  P  =  '^2Pj3qj,h  £  J\f+,qj  £  Qo,Pj  rational  ,^2pj  =  1} 

1=1  l=i 
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is  dense  in  (V(Q),p),  where  6qj  is  the  Dirac  measure  with  atom  at  qj. 

It  is  straight  forward  to  use  the  ideas  and  results  associated  with  this  theorem  to  develop  computationally 

OO 

efficient  schemes.  Given  Qd  =  \^j  QM  with  Qm  =  {Qj{}jLi  (a  “partition”  of  Q )  chosen  so  that  Qd  is  dense 

M= 1 

in  Q,  define 

M  k 

Vm(Q)  ={P  G  V(Q)  :  P  =  ^2pjSqM,qf  G  Qm,Pj  rational  .  ^ p?  =  1}. 


i= i 


j= i 


Then  we  find 

(i)  Vm(Q)  is  a  compact  subset  of  (V(Q),p), 

(\\)Vm(Q)  C  Pm+1(Q), 

(iii )“PM(Q)  — ►  P(Q)”  in  the  p  topology;  that  is,  elements  in  V(Q)  can  be  approximated  in  the  p  metric 
by  elements  of  VM  for  M  sufficiently  large. 

These  ideas  and  results  can  then  be  used  to  establish  a  type  of  “stability”  of  the  inverse  problem  (see 
(Banks  and  Bihari,  2001),  (Banks  and  Kunisch,  1989)).  We  first  define  a  series  of  approximate  problems 
consisting  of  minimizing 

n  2 

J(Pm)  =  ^2  |x(*b  Pm)  -  di 

i=  1 

over  Pm  G  Vm  (Q).  Then  we  have 

Theorem  5  Let  Q  be  a  compact  metric  space  and  assume  solutions  x(t\  P)  are  continuous  in  P  on  (V(Q),  p). 
Let  Qd  be  a  countable  dense  subset  of  Q  with  Qm  =  {qf1  }fLi  and  VM (Q)  as  above  so  that  (i)-(iii)  holds. 
Suppose  PM{dk )  is  the  set  of  minimizers  for  J{Pm )  over  Pm  G  Vm{Q)  corresponding  to  the  data  {dk}  and 
P*(d )  is  the  set  of  minimizers  over  P  G  V{Q)  corresponding  to  d,  where  dk ,  d  G  R"  are  the  observed  data  such 
that  dk  — >  d.  Then  dzst(P^(dfc),  P*(d))  — >  0  as  M  — >  oo  and  dk  —*  d.  Thus,  the  solutions  depend  continuously 
on  the  data  and  the  approximate  problems  are  method  stable  as  formulated  in  (Banks  and  Kunisch ,  1989). 

Of  course,  for  infinite  dimensional  state  systems  such  as  (1),  (3)  and  (4),  one  would  also  approximate 
the  solutions  x(t;PM )  by  finite  dimensional  approximate  solutions  xN (f;  Pm)  to  obtain  a  completely  finite 
dimensional  problem.  A  version  of  the  above  theorem  can  be  given  for  this  simultaneous  state/parameter 
approximation  using  the  approach  for  state/parameter  approximation  found  in  (Banks  and  Kunisch,  1989). 

The  “delta  measure”  approximations  given  above  are  essentially  zero-order  splines.  As  one  might  expect, 
higher  order  schemes  can  readily  be  developed.  An  example  of  linear  spline  schemes  has  been  developed  in 
(Banks  and  Pinter,  2005)  and  further  investigated  in  (Banks  and  Davis,  to  appear)  and  can  be  summarized 
in  the  following  theorem. 

Theorem  6  Let  IF  be  a  weakly  compact  subset  of  L2{Q),  Q  compact  and  let  Vr(Q)  =  {P  G  V(Q)  :  P'  = 
p,p  G  T }.  Then  Tr{Q)  is  compact  in  (V(Q),p).  Moreover,  if  we  define  {Ij1}  to  be  the  linear  splines  on  Q 

corresponding  to  the  partition  Qm,  where  [JQm  is  dense  Q,  and  define 

M 


VM  =  {pM  :pM  =  J2  rational  }. 

j 

Then  if 

VpM  =  {PM  G  V(Q)  :  Pm  =  J pM ,  pM  G  PM}, 
we  have  |^J  Pjc-m  is  dense  in  {Vr{Q),p). 

M 
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5  CONCLUDING  REMARKS 


The  framework  outlined  above,  while  most  useful,  is  not  complete.  Statistical  aspects  of  the  inverse  problems 
for  estimation  of  measures  as  discussed  above  are  still  under  investigation.  The  approximations  (delta  and 
splines)  lead  to  finite  dimensional  inverse  problems  for  which  standard  asymptotic  theory  (see  Chapter  12 
of  (Seber  and  Wild,  1989))  for  standard  errors  and  confidence  intervals  (using  the  sensitivity  functions  dis¬ 
cussed  above)  can  be  applied.  However,  an  analogous  asymptotic  theory  for  the  original  infinite  dimensional 
problems  involving  (5)  of  Section  4  has  yet  to  be  rigorously  developed. 
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